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PREFACE 


A recent  report,  based  on  a dissertation  prepared  by  the  first  author,  addressed  the 
problem  of  enhancing  the  accuracy  of  traditional  cubature  formulas  for  evaluating 
double  integrals  numerically  over  rectangles  by  the  addition  of  first-  and  mixed 
second-order  partial  derivative  correction  terms  evaluated  on  the  boundary  of  the 
domain  of  integration. 

The  present  report,  based  on  a paper  presented  at  the  Society  for  Industrial  and 
Applied  Mathematics  (SIAM)  1977  Fall  meeting,  31  October  to  2 November, 
Albuquerque,  New  Mexico,  generalizes  the  results  of  the  dissertation  to  multi- 
dimensional integrals  over  hyperrec tangles.  In  addition,  53  new  multidimensional 
quadrature  formulas  with  boundary  partial  derivative  correction  terms  are  given. 
Numerical  results  are  presented  for  double  and  triple  integrals. 
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MULTIDIMENSIONAL  QUADRATURE  FORMULAS  USING 
PARTIAL  DERIVATIVES 

1 2 
WAYNE  E.  HOOVER1  AND  J.  SUTHERLAND  FRAME 


I.  INTRODUCTION 

Traditional  methods  of  approximating  multidimensional  integrals  of  the  form 

1(f)  = f<xi>  * * * * xN,dxl  ' * • dxN 

over  the  hyperrec tangle 


r = n a.,  b.  | , ai*  bi  rea1’ 

i = 1 * 


employ  a weighted  sum  of  function  values 
N 

Q(f)  = £ w.f(x. ^ ...» 


The  w.  are  called  weights  and  the  (x.j,  . . . , *iN)  are  called  nodes.  The 
difference 

E(f)  = 1(f)  - Q(f) 

is  the  truncation  error  (or  error). 

Let  pk  = pk(xr  . . . ,xN)  be  a polynomial  of  degree  k in  N variables.  Then  the 

multidimensional  quadrature  formula  or  rule  Q(f)  is  said  to  be  of  order  k or  have 
degree  of  precision  k if  for  any  pk»  E(pk)  = 0,  but  E(pk+1)  * 0 for  at  least  one 

polynomial  Pk+j* 

Since  it  is  not  uncommon  for  numerical  procedures  to  make  use  of  partial 
derivatives,  e.g.»  in  optimization  techniques,  it  is  surprising  that  very  little  seems 
to  be  known  concerning  the  use  of  partial  derivatives  in  nonproduct  multi- 
dimensional quadrature  rules.  Stroud[21J  gives  only  one  cubature  formula,  C2:Z-1> 
due  to  Ionescue  [ll],  which  uses  partial  derivatives  of  the  integrand. 


^ Computer  Services  Directorate,  U.S.  Naval  Air  Test  Center,  Patuxent  River, 
Maryland  20670. 

^Department  of  Mathematics,  Michigan  State  University,  East  Lansing, 
Michigan  48824. 
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Therefore,  the  objective  of  this  investigation  is  to  construct  a number  of  new 
composite  multidimensional  quadrature  formulas  of  orders  3 and  5 using  first-  and 
mixed  second-order  partial  derivative  correction  terms  in  addition  to  function 
values  of  the  integrand.  Numerical  results  indicate  that  the  use  of  partial 
derivatives  evaluated  on  the  boundary  of  R increases  the  accuracy  and  efficiency 
of  composite  multidimensional  integration  formulas.  The  efficiency  of  the 
composite  formulations  may  be  explained  in  terms  of  the  following  three 
properties. 

First,  consider  an  m-point  rule  in  which  all  the  nodes  lie  in  the  interior  of  R. 
When  R is  partitioned  into  s subhyperrec tangles  or  cells  and  the  rule  is  applied  to 
each  cell,  the  total  number  of  nodes  is  ms  since  the  nodes  are  interior  to  each  cell. 
A more  efficient  procedure  is  to  employ  an  integration  rule  in  which  some  nodes 
coincide  with  the  boundary  of  the  domain  of  integration.  Then  when  the  domain  is 
subdivided,  these  nodes  are  included  in  more  than  one  cell.  Thus  the  total  number 
of  nodes  is  considerably  less  than  the  sum  of  their  numbers  in  each  cell.  This  we 
call  the  "inclusion  property." 

The  second  property  is  known  as  "persistence  of  form."  Briefly,  this  means  that 
for  some  functions  it  requires  approximately  the  same  if  not  less  computer  time  to 
evaluate  a partial  derivative  as  it  does  the  function  since  the  form  of  a partial 
derivative  follows  that  of  the  function.  Thus,  part  of  the  calculation  required  for 
the  function  may  be  reused  for  the  partial  derivative  evaluation.  Further  computer 
economy  can  also  be  achieved  if  the  partial  derivative  is  evaluated  at  the  same 
point  as  the  function  since  only  one  calculation  of  the  location  of  the  point  is 
required. 

Finally,  efficiency  results  from  applying  the  "alternating  sign  property." 
Essentially  this  means  that  in  the  composite  formulation  of  a rule,  the  weights 
assigned  to  the  partial  derivative  nodes,  equal  in  magnitude  but  opposite  in  sign, 
cancel  at  interior  points  and  consequently,  the  partials  need  be  evaluated  only  on 
the  boundary  of  R. 


In  Section  in,  the  method  of  undetermined  coefficients  in  conjunction  with  the 
inclusion  property,  the  persistence  of  form  property,  and  the  alternating  sign 
property  is  employed  to  construct  a number  of  new  partial  derivative  corrected 
multidimensional  numerical  integration  formulas. 


To  conclude  this  section,  observe  that  because  of  the  "boundary  effect,"  that 
is,  most  of  the  volume  of  a hyperrectangle  lies  near  the  boundary,  it  seems  natural 
to  construct  multidimensional  quadrature  rules  with  boundary  partial  derivative 
correction  terms.  Indeed,  for  the  hyperrec tangle  R with  edge  w^  = b^-a^,  the  set 
of  points  whose  distance  from  some  edge  is  less  than  or  equal  e/2  has  volume 


II  w,  - n (w.-ew.)  = I l-(l-e)N  I II  w 
k=l  k=l  1 ' k=l  K 
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which  approaches  the  volume  of  R as  the  number  of  dimensions,  N,  approaches 
infinity.  For  example,  the  volume  of  the  set  of  points  whose  distance  from  some 
edge  of  a 100-dimensional  hyperrectangle,  R,  is  less  than  or  equal  0.05  is  99.997% 
the  volume  of  R. 

H.  THE  ALTERNATING  SIGN  PROPERTY 

In  this  section,  several  observations  will  be  made  for  the  specific  case  N = 2. 

A study  of  the  Euler-Maclaurin  Summation  formula  for  a function  of  two 
variables  suggests  the  possibility  of  constructing  nonproduct  cubature  formulas 
involving  first-  and  perhaps  mixed  second-order  partial  derivative  correction  terms 
with  weights  of  equal  magnitude  but  alternate  signs  at  the  four  corners  or  at  the 
midpoints  of  the  sides  of  the  rectangular  domain  of  integration, 
R = [ -hj,  hj  ] x [ -h.,,  h^],  so  that  when  the  rule  is  compounded  or  repeated,  the 
weights  cancel  except  on  the  boundary.  As  previously  noted,  this  is  called  the 
alternating  sign  property.  See  Figure  1. 


Figure  1 

Alternating  Sign  Property  for  N = 1 and  2 
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This  investigation  will  consider  the  six  "cubature  elements"  listed  in  Table  I. 

Table  I 

Cubature  Elements 


Centroid 

Value 

Vertex 

Sum 


Midpoint 

Sum 


Vertex 

Derivative 

Correction 


Midpoint 

Derivative 

Correction 


Vertex 

Mixed 

Derivative 

Correction 


Cubature  Element 


Diagram 


Qc(f)  = -Ihjh^f (o,  o) 


Qv(f)  = [f(+,+)*  + f(-,+)  + f(-,-)  + f(+,-)] 

Qm(f)  = [«♦*•  + + f(-,o)  + f(o,-)j 


Qv  (f»  - 4h?h2  [fx(*’*'  - * 'x1*'-']  .» 1+ 

* 4h,h*  [fy(»,*l  * fy<-,*)  - - 


Q;ifl  - 4hihz  [yso,  - f„(-,oi] 

* 4hjh^  [fy(0,.|  - fy(o,-)] 
Q;tt)  . 4hjh^  - fxyW) 


* 'xy'-'-1  - V'->] 


-Q 

□ 


*f(+,+)  = f (+h  j , +h2) 


The  selection  of  the  derivative  correction  cubature  elements  O'  , O'  , and  Q" 

v m v 

is  based  on  the  following  considerations.  For  a and  6 nonnegative  integers,  define 
the  partial  derivative  correction  terms,  T.,  listed  in  Table  II. 
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Table  H 

Partial  Derivative  Correction  Terms 


Now  suppose  flxpX^)  can  be  expanded  in  a Taylor  series  about  the  point  (0,0) 

as  far  as  may  be  required.  If  the  series  converges  on  R,  then  the  correction  terms 

assume  zero  or  nonzero  values  as  indicated  in  Table  ID.  Examination  of  Table  III 

shows  the  odd/even  constraints  which  must  be  placed  upon  a and  8 in  order  to  avoid 

nonzero  correction  terms,  T.,  the  components  of  the  partial  derivative  correction 

elements  Q'  Q'  and  Q"  . These  cubature  elements  are  then  candidates  for 
v m v 

inclusion  in  cubature  formulas.  In  Section  VI,  the  n-dimensional  generalizations  of 
the  elements  are  used  as  building  blocks  to  construct  53  new  multidimensional 
quadrature  rules  with  partial  derivative  correction  terms.  The  results  are 
summarized  in  Table  VII. 


Name 


Correction  Term 


T5  [«+.+)-«-.+)  + «-r)  - f(+r) 


Diagram 
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f 


l 


Table  in 

Values  of  the  Correction  Terms 


Correction 

Term 


a Even 
8 Even 


T 


1 


0 


a Odd 
8 Even 


Nonzero 


a Even 
8 Odd 


0 


a Odd 
8 Odd 


0 


0 Nonzero 


Nonzero  0 


0 


0 


0 Nonzero  0 


0 0 Nonzero 


ID.  DERIVATION  OF  MINTOV 

Denote  the  center  of  the  hyperrectangle 


R 


■ ,?i  M 


N 

by  c = (0,.  . . ,0),  the  2 vertices  by  v = (vj,  . • . ,v^)  where  v.  = h.  or  -lv,  the  2N 

midpoints  of  the  bounding  (N- l)-dimensional  hyperrectangles^  or  ^sides'' "*of  R by 

m = (m , , . . ,m. ,)  where  m.  = h.  or  -h.  and  m.  = 0 for  i i i,  and  the  volume  by 
I N j j j i J 7 

N 

h = n 2h.. 


j = l 


J 


For  p = m or  v,  define  the  sign  functionals 


^ + 1 ot! 

ojk(p>  = o (p)  ak(p) 


Oj(p)  = J -1  if  Pj  = ~hj 
otherwise 


(1) 


6 


“ — 


4 


and  the  first-  and  second-order  partial  derivative  correction  terms 


E»jf(P)  = J o (p)(J)'  f(p) 

(2) 

V (v)  = l « (v)(/;  (J) • f(v) 
where  the  sums  are  over  the  indicated  points  and 


denotes  the  a-th  partial  derivative  of  f(x)  with  respect  to  the  j-th  variable. 
Figures  2 and  3 illustrate  two  of  the  partial  derivative  correction  terms  for  the 
4-cube.  Careful  inspection  will  reveal  the  sign  arrangements  of  D^lv)  and  D1?f(v) 
for  the  2-  and  3-  cubes. 
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Figure  3 

Sign  Arrangement  for  Dj.,f(v) 

Now,  since  1(f)  = / f and  the  generalized  cubature  elements 

' R 


Qc(f>  - 

hf(c) 

Qp(f)  = 

h p tipi 

N 

Q’Jf)  = 

h E h.D.f(p) 

P 

j = l ] J 

Qy(f)  = 

» j<khAD: 

f (v) 


(3) 


vanish  for  functions  which  are  odd  in  any  variable,  we  may  approximate  1(f)  by  the 
linear  combination 


1(f)  - Q(f)  = X]Qc  + X2Qv  + X3Q'v  + X4Q'^ 


(4) 


• « • , , ,2422642  . 2 2 2 . 

which  is  exact  for  the  even  functions  1,  Xj,  Xj,  x^,  Xj,  x^,  and  xjx2x3-  By 

symmetry  then,  Q(f)  will  also  be  exact  for  all  polynomials  of  degree  at  most  5. 
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THEOREM.  If  f(x)  has  continuous  partial  derivatives  of  the  first  six  orders  on 
Ur  i 

R = II  I -h.,h.  . then 

i=i  L J 

1(f)  = |^8Qc  + (7Qv  - Q'y  - Q^/3)/2N  j /15  + E(f)  (5) 

is  a multidimensional  quadrature  formula  with  degree  of  precision  5.  The  truncation 
error,  E(f),  is  bounded  by 

| E(f)  | <_h  Je6  + 35E42  + 280E222]  /9450  (6) 


where 


4a8 . . .y 
jk.  . . L max 


<f>]  'w 


N a a 
E h.  M. 

3 = 1 3 3 


: = E hah®  MaS 

a6  j>k=1  3 k 3k 


v = l ua.S.  Ma87 
ot 87  j<k<L  j k 1 jkL  • 

PROOF.  By  applying  Taylor's  theorem  for  N-variables  to  (4)  and  equating 
coefficients  of  similar  terms,  one  obtains  the  linear  system 


0 O X 


1 _N-3  1 ,N-1 


X 
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having  the  unique  solution 

(Aj,  X2*  A^i  A^)  = 


(±  JL d d_\ 

\ 15  ’ 2N 1 5 ’ 2N15  2N45  / 


(9) 


The  bound  (6)  on  the  truncation  error  is  a consequence  of  Taylor's  theorem  and  is  a 
straightforward  calculation.  Observe  that  the  last  term  in  the  error  bound  appears 
only  for  N > 2. 

The  name  MINTOV  (for  Multiple  INTegration,  Order  5)  is  given  to  (5). 

IV.  COMPOSITE  FORMULATION  OF  MINTOV 

In  order  to  obtain  the  composite  formulation  of  (5)  for  an  arbitrary 
hyperrectangle 


partition  each  interval 


■ l M' 

[v>i] 


into  n.  subintervals  each  of  length  h.  = w./n.  where 


] J J 


w.  = b.  - a..  Denote  the  volume  of  each  cell  thus  obtained  by 
J J J 

N 

h = n h. 


j=l 


J 


and  the  volume  of  R by 
N 

w = II  w 


j = l 


I 


To  condense  notation,  define 

U0) 


(aj+hj  (ij-0),  • • • , aN+hN(iN-0)) 


C(x.,0) 


C(x.) 


(aj+hjUj-0),  . . . , x.,  a.  + 1+h.  + 1(L  + 1-0) 


aN  + hN(,N'0)) 


(10) 


= (aj+ijhj,  . . . , x.,  aj  + 1 + ‘j+1hj+i»  • ‘ • * aN+,NhN^ 
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S(x.,xk)  = Uj+ijhj,  ....  x.,  aj+1  +ij+1hj+1,  . . • xR, 


aN+iNhN) 


and 


N a 

E • E 


°N 


n2  nl 


a=l  ia=l 


*nb1 


*2  1 'r1 


Furthermore, 


n 


N a 

Ac  = h E • z f[UH)] 


a=l  ,a=l 


n 


N "a 
Av  = h E • E 7f[U0)] 


a=l  ia=0 


(11) 


m 


N 

h E 

j = l 


N 

E • 

a = l 


1 rf[c,ar,h.,«] 

‘a-  ‘a  l, 


A' 


N 

h E 

j=l 


N n„ 


E • E 


a = l ia=0 


7h.  {/)'.  [f(C(b.))-f(5(a.))] 


A ' 
m 


N 

h E 

j=l 


N 

E 

a = l 


a 

E h (J)  ' [f(C(b  ,H))  - f(C(a.,  %))] 
i =1  J J J J 


N 

E 


A;  = h E 

j<k  a = l ia=0 


no 

E 7hjhk^k^>j  Cf<C(aj,ak))  - £(C(a.,bk)) 


j.k 


- f(C(b.,ak))  + f(5(b.,bk))] 
Ooserve  that  the  weight  7 is  to  be  assigned  when  the  node  is  common  to 7 cells. 
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I 


COROLLARY.  If  f(x)  has  continuous  partial  derivatives  of  the  first  six  orders  on 
the  hyperrec  tangle  R,  then 


1 

f(xj,  ...  jX^ldxj  ...  dx^j 

"N  “1 

= [ 8Ac  + (7Ay  - A'y/2  - A;/12)/2N  ) /15  + E(f) 


(12) 


is  a composite  multidimensional  quadrature  formula  with  first-  and  second-order 
partial  derivative  correction  terms  having  degree  of  precision  5.  The  truncation 
error  is  bounded  by 

| E(f)  | <_  w |e6  + 35E42  + 280E222  j /604  800.  (13) 

PROOF.  Let  v = (vj,  . . . ,vN),  v.  = a.  or  b.,  represent  a vertex  of  R,  c = 
(Cj,  . . . , cN),  Cj  = Vi(a.  + b.),  the  cell  centroid,  and  m = bj,  . . . , mN>,  nr  = a.  or 
b.  and  nr  = c.  for  i^  j,  a midpoint  of  a side.  For  p = m or  v define  the  sign 
functionals 


Oj(p) 


{ 


-1  if  Pj  = a, 

+ 1 otherwise 


o.k(p)  = a^(p)ak(p). 


(14) 


Then  by  a linear  change  of  variables,  using  the  notation  of  (2),  and  replacing  h and 
Ir  in  (3)  by  w and  w.,  respectively,  one  obtains 


r bN  rb i 

/ • • • / f(tj,  • • • ’^n^I  ■ * ■ dtN 

J aN  al 


/ •••/ 


‘N 


f( 


W1X1 

2h, 


WNXN.  w . 

•’-Zh^h  dxl-**dxN 


(15) 


I 


12 
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f 


k 


I 


_8_ 

15 


Qc  + N 
C 2 15 


— O' 

2N+115  V 


1 


2N+245 


Q"  + E(f). 
v 


The  proof  follows  by  applying  (15)  to  each  subhyperrec tangle  of  R.  Here  again,  the 
last  term  in  (13)  appears  only  in  the  case  N > 2. 

The  name  MINTOV  is  also  given  to  (12)  since  it  is  the  composite  formulation  of 
(5).  With  appropriate  interpretation  of  the  centroid  and  corner  nodes  and  the  partial 
derivative  correction  terms,  MINTOV  is  a nonproduct  N-dimensional  generalization 
of  the  composite  Simpson's  formula  with  end  corrections  (Lanczos,  [12]): 


f (x)dx 


8h  “ 

15  L 
3 i=l 


f(a+h(i-  H))  + 


n 

E ' f(a+ih) 
i=0 


H - *■<«>]  * sssioo  h6,l6’|^,• 


(16) 


The  prime  on  the  summation  signifies  that  the  weight  7 is  to  be  assigned  in  case 
the  node  a+ih  is  common  to  7 subintervals.  Also,  n=(b-a)/h,  and  £ is  some  point  in 
£ a,b  ] . The  number  of  function  evaluations  is  2n  + 3. 

It  may  also  be  stated  that  MINTOV  is  composite  Ewing's  formula  [ 5 J with  partial 
derivative  correction  terms.  In  Section  VI,  an  N-dimensional  composite  formulation 
of  Ewing's  formula  is  given  as  D0503. 

V.  NUMBER  OF  FUNCTION  EVALUATIONS 

In  this  section  a partial  derivative  evaluation  is  counted  the  same  as  a function 
evaluation.  For  many  integrands  it  may  be  necessary  to  weight  the  partial 
derivative  evaluations.  However,  because  of  the  persistence  of  form  property 
previously  noted,  the  present  enumeration  technique  will  be  used. 

Indeed,  for  functions  such  as  ln(xyz)  the  first-order  partials  require  only  about 
40%  of  the  time  to  evaluate  the  given  function,  whereas  functions  similar  to 

cos(x)cos(y)cos(z),  exp(-xyz),  (1+x+y+z)-4,  and  even  (l+w)sin(x)sin(y)sin(z)e”W/ 

2 2 2 2 

(xyz),  w = x +y  +z  , require  not  more  than  5%  additional  time  to  evaluate  the 
first-  and  mixed  second-order  partials  than  the  original  functions. 
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Now  the  numbers  of  function  evaluations,  p , required  by  the  functionals  in 
(11)  are 


P(AC) 


N 

' J, 


N 

p(A  ) = II  (n.  + l) 

i=l 


P(AJ 


N N 

Z (n.  + l)  fl  n. 


3=1 


3 


i=l 

i/t  j 


p(A'v)  = 


N N 

Z Z n (n.  + l) 
3=1  i=l  1 


(17) 


N N 


p(A'  ) = 2 Z 

m 


n 


n. 


3=1  »=1 

^3 


P(A" 


N 

4 £ I!  (n.  + l). 
j<k  i=l 


Thus,  for  d-dimensional  MINTOV,  if  n^  = (b.  - a.)h.  = n for  all  i,  the  number  of 
function  evaluations  is 

p = nd  + (n+l)d  + 2d(n+l)d_1  + 2d(d-l)(n+l)d'2 


= 2n 


+ V [^2(d-i)2  + l]  ( f ) n\ 


(18) 


Hereafter,  we  refer  to  n as  the  number  of  partitions  of  R.  Note  that  R is 
subdivided  into  n subregions. 
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Let  Lyness,  Gauss,  and  Boole  represent  the  composite  formulations  of  Cn:5-5 

as  listed  in  Stroud  [ 2l]  and  which  is  due  to  Mustard,  Lyness,  and  Blatt  [ 16^  , the 
composite  product  Gauss,  and  the  fifth-order  composite  product  Newton-Cotes 
formulas,  respectively.  The  number  of  function  evaluations  for  each  rule  is  given  in 
Table  IV.  For  example,  for  eight  partitions  in  4-space,  MINTOV,  Lyness,  Gauss,  and 
Boole  require  18  433,  43  425,  331  776,  and  1 185  921  function  evaluations, 
respectively. 

Table  IV 

Number  of  Function  Evaluations  for  Several  Fifth-Order  Formulas 


Formula 

P 

MINTOV 

nd  + (n+l)d'Z  (n+1)2  + 2d(n+d) 

LYNESS 

(n+1 )d  + (2d+l )nd 

GAUSS 

(3n)d 

BOOLE 

(4n+l)d 

It  can  be  shown  that  for  n >d,  p (MINTOV)  < p (Lyness);  equality  holds  only  for 
n=d=2. 

Of  special  interest  are  the  cases  d = 2 and  3.  In  the  case  d = 2,  observe  that  the 
composite  formulation  of  the  9-point,  degree  5 Lyness  cubature  formula  requires 

6n2  + 2n  + 1 function  evaluations,  fe.  The  Radon  fl8]  , Albrecht,  Collatz  [ l} 

7-point,  degree  5 composite  formula  requires  7n2  fe.  For  brevity,  we  call  it  Radon's 
formula. 

As  Table  IV  shows,  the  9-point,  degree  5 composite  product  Gauss  cubature 
formula  requires  9n2  fe,  and  the  25-point,  degree  5 composite  product  Boole's 
(Newton-Cotes')  rule  requires  16nZ  + 8n  + 1 fe. 

Tanimoto's  [ 22 } fifth-order  derivative  corrected  Simpson's  rule  requires 

4n2  + 12n  + 9 fe,  whereas  the  fifth-order  MINTOV  and  DH5G5S  (see  Table  VII) 
2 2 

rules  require  only  2n  + 6n  + 9 fe  and  2n  + lOn  + 9 fe,  respectively.  The  composite 
product  formulation  of  the  third  order  Simpson's  rule  requires  4n  + 4n  + 1 fe. 

Table  V shows  the  number  of  function  evaluations  these  cubature  formulas 
require  for  various  subdivisions.  Similar  results  for  the  case  d = 3 are  presented  in 
Table  VI.  As  previously  noted,  a partial  derivative  evaluation  is  counted  as  one 
function  evaluation. 
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Table  V 

Number  of  Function  Evaluations  for  Various  Cubature  Formulas 

Compounded  n^  Times.  These  are  Fifth-Order  Formulas  Except  for 
Simpson's  Rule  which  is  Third-Order 


■ 

MINTOV* 

Zn2+6n+9 

SIMPSON 

4n^+4n+l 

TANIMOTO* 

4n^+12n+9 

LYNESS 

6n^+2n+l 

RADON 

7„2 

GAUSS 

9n2 

BOOLE 

I6n^+8n+l 

1 

17 

21 

9 

25 

mm, 

7 

25 

2 

29 

37 

25 

49 

fjfl 

28 

IB 

81 

4 

65 

81 

81 

121 

105 

112 

144 

289 

8 

185 

217 

289 

361 

401 

448 

576 

1 089 

16 

617 

681 

1 089 

1 225 

1 569 

1 792 

2 304 

4 225 

32 

2 249 

2 377 

4 225 

4 489 

6 209 

7 168 

9 216 

16  641 

64 

8 585 

8 841 

16  641 

17  161 

24  705 

28  672 

36  864 

66  049 

100 

20  609 

21  009 



40  401 

41  209 

60  501 

70  000 

90  000 

160  801 

•Partial  Derivative  Corrected  Formulas 


Table  VI 

Number  of  Function  Evaluations  for  Various  3-Dimensional  Quadrature 

3 

Formulas  Repeated  n Times 


r 

MINTOV* 

2n349n2*27n*l9 

DH5G5S* 

2n3*15n2427ifl9 

SH9G5R* 

Sn3*l8n2427n»l9 

LYNESS 

8n5.)nZ*Jn«l 

SIMPSON 

8n3* I2n2*6n* I 

GAUSS 

27n3 

BOOLE 

64n3*48n2* 12n*  1 

B 

57 

63 

69 

15 

27 

27 

64 

1 

125 

149 

185 

83 

125 

216 

729 

B 

399 

495 

735 

373 

729 

I 728 

4 913 

B 

1 835 

2 219 

3 947 

4 313 

4 913 

13  824 

35  937 

16 

10  947 

12  483 

25  539 

33  585 

35  937 

110  592 

274  625 

32 

75  635 

81  779 

183  155 

265  313 

274  62$ 

884  736 

2 146  689 

64 

562  899 

587  475 

1 386  195 

2 109  633 

2 146  689 

7 077  888 

16  974  593 

100 

2 092  719 

2 152  719 

5 182  719 

8 030  301 

8 120  601 

27  000  000 

64  481  201 

‘Partial  Derivative  Corrected  Formulas  (cf.  Table  VID 
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VI.  CONSTRUCTION  OF  53  NEW  MULTIDIMENSIONAL  QUADRATURE  RULES 
WITH  PARTIAL  DERIVATIVE  CORRECTION  TERMS  AND  ERROR 
ESTIMATES 

As  in  the  case  of  MINTOV,  the  method  of  undetermined  coefficients  is  used  to 
construct  a number  of  new  formulas.  The  results  are  compiled  in  Table  VII;  they  are 
n-dimensional  generalizations  of  the  cubature  formulas  given  in[10].  Observe 
that  the  absolute  values  of  the  entries  in  the  Error  columns  are  to  be  used  for  the 
error  estimates.  The  minus  signs  are  a consequence  of  applying  the  method  of 
undertermined  coefficients  and  are  included  for  comparison  purposes.  Also,  note 
that  some  of  the  terms  become  zero  for  n = 3,  4,  and  7. 

Entry  26,  DC5C5,  is  MINTOV  as  given  in  (12)  and  (13).  For  ease  of  reference, 
the  2-dimensional  formulation  will  be  stated  explicitly.  To  this  end  let  the 
rectangle  R = [a,b]  x [ c,d  ] be  partitioned  into  nm  cells  each  of  size 
hk  = (b-a)(d-c)/nm.  Then 


/:/ 


f (x,y)dxdy  = E E f[a+h(i-H),  c+k(j-H)] 


m n 


15 


7hk 


j=l  i=l 


m n 


E'  E’  f(a+ih,  c+jk) 
j=0  i=0 


,2,  m 

h k r i 

120  jfo 


hk 


2 n 


120  . 


r 


i*0 


I f (a+ih,d)  - 

Ly 


c+jk)  - fx(a,c+jk) 


(19) 


f^(a+ih,c) 


h\2  r 


720 
+ E(f) 


f (a,c)  - f (b,c)  + f (b,d)  - f (a,d) 
xy  xy  xy  xy 


p (MINTOV)  = 2 (nm)  + 3(n+m)  + 9. 


(h6Mj  + kfeM2>  + 35(h4k2M42  + h2k4M24)  (20) 


(21) 


The  selection  of  the  formula  names  assigned  to  the  formulas  listed  in  Table  VII 
was  based  on  the  following  considerations.  The  first  two  symbols  were  chosen 
somewhat  arbitrarily,  whereas  a digit  was  selected  for  the  third  symbol  to 
represent  the  number  of  function  evaluations  required  for  the  holistic  cubature 
rule,  that  is,  for  d = 2 and  n^  = n^  = 1. 
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Table  VII  (Cont'd) 
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The  fourth  symbol  represents  the  number  of  partial  derivative  evaluations 
required  for  the  holistic  cubature  rule.  The  fifth  digit  is  the  degree  of  precision  of 
the  formula.  The  presence  of  a sixth  symbol  signifies  that  the  method  of 
undetermined  coefficients  led  to  more  equations  than  unknowns.  The  letters  A and 
B are  used  to  indicate  two  successful  combinations,  whereas  an  R,  S,  or  T indicates 
at  least  one  unsuccessful  combination. 

For  completeness,  formulas  1 and  7,  the  composite  formulations  of  the 
midpoint  [93  and  trapezoidal  rules,  respectively,  have  been  included.  Formula  51 
is  equivalent  to  the  composite  Simpson's  rule  for  N=2.  The  holistic  representation 
(i.e.,  na=  1 for  all  a)  for  formulas  13,  19,  30,  and  40  were  investigated  by  Squire 
1 20 J , Ewing  [5],  Tyler  [23],  and  Miller  [15],  respectively.  Formula  30  was 
apparently  discovered  independently  by  Bickley  [2]  and  Tyler  [23].  We  will  follow 
Stroud  [ 2 1 3 and  call  it  Tyler's  rule.  Thus,  except  for  formulas  1,  7,  13,  19,  30,  40, 
and  51,  the  remaining  53  multidimensional  quadrature  formulas  are  new. 

In  the  two-dimensional  case,  it  is  interesting  to  compare  DF543S  with 
Simpson's  cubature  rule,  S0903S,  because  both  are  third-order  formulas  and  both 
have  the  same  error  bounds;  however,  the  former  requires  approximately  half  as 
many  function  evaluations  as  Simpson's  rule.  Specifically,  the  formulas  require 
2(njn^)  + (n^+n^)  + 5 and  4(njn^)  + 2(nj+n^)  + 1 function  evaluations, 
respectively. 

DF543S  is  a combined  trapezoidal-midpoint  or  Ewing  rule  [ 5]  with  boundary 
correction  terms.  This  cubature  rule  requires  mixed  second-order  partial 
derivatives  evaluated  only  at  the  corners  of  the  rectangular  domain  of  integration 
R. 

Two  formulas,  XF543S  and  OF893S,  which  share  this  property  with  DF543S 
also  have  the  same  error  bound  as  Simpson's  rule.  There  may  be  applications  where 
the  cubature  rule  DF543S  should  be  considered  as  a viable  alternative  to  Simpson's 
rule. 


Finally,  we  note  that  numerical  results  for  double  integrals  indicate  that 
EM143  is  competitive  with  Simpson's  rule  since  it  requires  only  (n^n^)  + 2(nj+n2> 
function  evaluations  as  compared  with  ^n^)  + 2(n1+n2>  + 1 for  Simpson. 
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Numerical  results  are  presented  in  Table  VIE.  As  before,  p represents  the 
number  of  function  (and/or  partial  derivative)  evaluations. 

Table  VIII 


— 

yoyo  1 + x 


dxdy 

j 2 2 

1 + x y 


= 0.915  965  594  177 


h = k - TO 
(n=n>  = 10) 


MIDPOINT 

E01P 

EM  143 

ET183 

EX183S 

EC1C3S 

ES1C3S 


TRAPE- 

ZOIDAL 

TP4P1 

TMr43 

TT433 

TX483S 

TC4C3S 

TS4C3S 


SQUIRE 

M(J4P1 

MM443 

MT483 

MX483S 

MC4C3S 

MS4C3S 


EWING 

D05P3 

DF543S 

DM  543  A 

DM543B 

DT583A 

DT583B 

DX585 

MINTOV 

DS5C5 

DH5G5S 


Error 

-9 

. 52-4^ 

-1 

.11-6 

-1 

.15-6 

-1 

.11-6 

-1 

.15-6 

-1 

.11-6 

1 

.90-3 

1 

.18-6 

1 

.27-6 

1 

.18-6 

1 

.27-6 

1 

.24-6 

3 

.76-4 

1 

.01-7 

1 

.22-7 

1 

.01-7 

1 

.22-7 

9 

.39-8 

-3 

.44-7 

-3, 

.44-7 

-8 

.53-7 

-3 

.88-8 

1 

.93-7 

-2 

.20-8 

-3 

.88-8 

-2 

.20-8 

-I, 

.64-8 

4, 

.31-10 

Error 

Order 

-2.38-4 

1 

-6.97-8 

3 

-7.04-8 

3 

-6.97-8 

3 

-7.04-8 

3 

-6.98-8 

3 

4.76-4 

1 

7.84-8 

3 

7.97-8 

3 

7.84-8 

3 

7.97-8 

3 

7.93-8 

3 

1.19-4 

1 

5.32-9 

3 

5.64-9 

3 

5.32-9 

3 

5.64-9 

3 

5.21-9 

3 

-2.04-8 

3 

-2.04-8 

3 

-5.33-8 

3 

-6.01-10 

3 

1.30-8 

3 

-3.39-10 

3 

-6.01-10 

5 

-3.39-10 

5 

-2.52-10 

5 

8.61-12 

5 
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Table  VIII  (Cont'd) 


IB 

B 

No. 

Formula 

P 

Error 

p 

Error 

Order 

30 

TYLER 

X0503 

85 

-3.02-7 

320 

-1.97-8 

3 

31 

XF543S 

89 

-3.02-7 

324 

-1.97-8 

3 

32 

XM543T 

105 

2.04-8 

360 

3.14-10 

3 

33 

XT583A 

109 

-4.43-7 

364 

-2.81-8 

3 

34 

XT583B 

109 

3.72-8 

364 

5.75-10 

3 

35 

XX585 

109 

2.04-8 

364 

3.14-10 

5 

36 

XC5C5 

113 

3.72-8 

368 

5.75-10 

5 

37 

XS5C5 

129 

1.34-8 

404 

2.06-10 

5 

39 

XH5G5S 

133 

7.31-10 

408 

9.76-12 

5 

40 

MILLER 

00803 

n 

-2.60-7 

341 

-1.90-8 

3 

41 

OF843S 

-2.60-7 

345 

-1.90-8 

3 

42 

OM843A 

mi 

2.21-7 

381 

1.34-8 

3 

43 

OM843B 

2.88-8 

381 

4.45-10 

3 

45 

OT883B 

mi 

4.56-8 

385 

7.06-10 

3 

46 

0X885 

2.88-8 

385 

4.45-10 

5 

47 

OC8C5 

124 

4.56-8 

389 

7.06-10 

5 

48 

OS8C5 

140 

1.76-8 

425 

2.71-10 

5 

50 

OH8G5S 

144 

7.74-10 

429 

9.92-12 

5 

51 

SIMPSON 

S09P3S 

121 

-3.16-7 

441 

-1.99-8 

3 

52 

SM945 

141 

6.27-9 

481 

9.65-11 

5 

53 

ST985 

145 

-1.07-8 

485 

-1.65-10 

5 

55 

SX985S 

145 

6.31-10 

485 

9.37-12 

5 

57 

SC9C5S 

149 

5.46-10 

489 

9.04-12 

5 

59 

SS9C5S 

165 

6.03-10 

525 

9.26-12 

5 

TANIMOTO 

169 

5.91-10 

PI 

9.22-12 

5 

LYNESS 

161 

5.66-9 

8.70-11 

5 

RADON 

175 

-1.84-9 

700 

-2.81-11 

5 

GAUSS 

225 

1.78-10 

900 

2.83-12 

5 

BOOLE 

441 

-1.85-10 

1681 

-2.77.12 

5 
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These  examples  were  run  on  the  U.S.  NAVAIRTESTCEN's  Real-Time 
Telemetry  Processing  System  Xerox  Sigma  9 computer  using  the  CPR  version  E00 
operating  system  and  FORTRAN  IV  in  double  precision  with  15+  significant  decimal 
digits. 

Example  2 (Stroud,  ^21^  ) 


/I  /•  1 j 

/.,  \ "" y 


3+x+y  dxdy  = (1  - 18  VT  + 25  5 ) 


= 6.859  942  640  334  65. 


(24) 


Results  for  h-k-^  are  presented  in  Table  IX.  For  this  example,  the  computed 
error  bounds  provide  reasonably  close  estimates  to  the  actual  errors. 


Figure  5 

Graph  of  z = V 3+x+y  on  f-l,lj  2 
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Table  IX 


3+x+y  dxdy  = 6.859  942  640  334  65 


h 

= k = ^ (n  = m 

= 6) 

* 

No. 

Formula 

P 

Error 

Error 

Bound 

Order 

i 

MIDPOINT 

E0101 

36 

-2.10-3+ 

9.26-3 

1 

2 

EM143 

60 

1.44-6 

1.93-4 

3 

3 

ET183 

64 

2.38-5 

5.14-4 

3 

4 

EX183S 

64 

5.21-6 

1.13-4 

3 

5 

EC1C3S 

68 

4.96-6 

1.13-4 

3 

6 

ES1C3S 

88 

5.17-6 

1.13-4 

3 

7 

TRAPE- 

ZOIDAL 

T0401 

49 

4.23-3 

1.85-2 

1 

8 

TM443 

73 

2.73-5 

7.72-4 

3 

9 

TT483 

77 

-2.11-5 

4.50-4 

3 

10 

TX483S 

77 

-6.47-6 

1.29-4 

3 

11 

TC4C3S 

81 

-5.96-6 

1.29-4 

3 

12 

TS4C3S 

101 

-6.13-6 

1.29-4 

3 

13 

SQUIRE 

M0401 

84 

1.05-3 

4.63-3 

1 

14 

MM443 

108 

-4.03-6 

8.84-5 

3 

15 

MT483 

112 

-1.52-5 

3.30-4 

3 

16 

MX483S 

112 

-2.51-7 

8.04-6 

3 

17 

MC4C3S 

116 

-1.22-7 

8.04-6 

3 

18 

MS4C3S 

136 

-2.94-7 

8.04-6 

3 

19 

EWING 

D0503 

85 

8.87-6 

1.93-4 

3 

20 

DF543S 

89 

1.32-6 

3.22-5 

3 

21 

DM543A 

109 

3.91-6 

8.57-5 

3 

22 

DM543B 

109 

1.18-5 

2.57-4 

3 

23 

DT583A 

113 

-1.11-6 

2.14-5 

3 

24 

DT583B 

113 

2.88-6 

6.43-5 

3 

25 

DX585 

113 

-2.41-7 

1.67-5 

3 

26 

MINTOV 

117 

-1.38-7 

9.65-6 

5 

27 

DS5C5 

137 

-1.04-7 

7.30-6 

5 

29 

DH5G5S 

141 

-6.98-10 

2.68-7 

5 

•See  Table  vn 

'-2. 10-3  = 2.10  x 10 
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Table  IX  (Cont'd) 


h = 

= k = i (n  = m 

= 6) 

No. 

Formula 

P 

Error 

Error 

Bound 

Orde 

30 

TYLER 

X0503 

120 

-2.21-6 

1.13-4 

3 

31 

XF543S 

124 

1.57-6 

3.22-5 

3 

32 

XM543T 

144 

-3 . 66-6 

8.04-5 

3 

33 

XT583A 

148 

2.13-6 

4.55-5 

3 

34 

XT583B 

148 

-1.26-5 

2.73-4 

3 

35 

XX585 

148 

1.13-7 

8.47-6 

5 

36 

XC5C5 

152 

2.16-7 

1.55-5 

5 

37 

XS5C5 

172 

7.05-8 

5.54-6 

5 

39 

XH5G5S 

176 

-6 . 66-9 

2.68-7 

5 

40 

MILLER 

00803 

133 

-1.33-5 

3.54-4 

3 

41 

OF843S 

137 

1.83-6 

3.22-5 

3 

42 

OM843A 

157 

-9.42-7 

2.14-5 

3 

43 

OM843B 

157 

-5.88-6 

1.29-4 

3 

45 

OT883B 

161 

-1.48-5 

3.22-4 

3 

46 

0X885 

161 

1.64-7 

1.20-5 

5 

47 

OC8C5 

165 

2.67-7 

1.90-5 

5 

48 

OS8C5 

185 

9.53-8 

7.30-6 

5 

50 

OH8G5S 

189 

-7.51-9 

2.68-7 

5 

51 

SIMPSON 

S0903S 

169 

1.49-6 

3.22-5 

3 

52 

SM945 

193 

2.90-8 

2.61-6 

5 

53 

ST985 

197 

-7.04-8 

4.96-6 

5 

55 

SX985S 

197 

-4.67-9 

2.68-7 

5 

57 

SC9C5S 

201 

-2.97-9 

2.68-7 

5 

59 

SS9C5S 

221 

-4.10-9 

2.68-7 

5 

TANIMOTO 

225 

-3.88.9 

_ 

5 

LYNESS 

229 

3.28-8 

- 

5 

RADON 

252 

-1.12-8 

- 

5 

GAUSS 

324 

-1.16-9 

- 

5 

BOOLE 

625 

1.21-9 

- 

5 

2' 
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Example  3 
ft  ft 

/:/:/ 


sin(x)sin(y)sin(z)  e*W  dxdydz 

xyz 


(25) 


= 1.531  670  226  93, 

where  w2  = x2+y2+z2.  The  first-order  partial  derivative  with  respect  to  x is 


r 2 sin(x)  1 sin(y)sin(z)e~w 

(x,y,z)  = (l+w)cos(x)  - (1+w+x  ) x J xyz 

and  the  mixed  second-order  partial  with  respect  to  x and  y is 

r 2 ,2  2.  2 2 

I w+w  +w(x  +y  ) + x y 

f (x,y,z)  = |(l+w)cos(x)cos(y)  + sinixjsimyi 

1 sin(z)e  W 

J xyz 


(26) 


2 2 

sin(x)cos(y)  - cos(x)sin(y) 

x y 


(27) 


The  results  of  applying  each  of  the  formulas  of  Table  VII  to  the  integral  in  (24) 
with  n.  = 8(h.  = */16>  are  presented  in  Table  X.  Note  that  for  triple  integrals,  the 
following  formulas  listed  in  Table  VII  coincide:  13=14=15;  16-17;  31-41;  and 
48=50=59. 
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vm.  CONCLUSIONS 

We  have  investigated  the  problem  of  enhancing  the  accuracy  of  conventional 
formulas  for  evaluating  multiple  integrals  numerically  over  n-dimensional 
rectangles  by  the  addition  of  first-  and/or  mixed  second-order  partial  derivative 
correction  terms  evaluated  on  the  boundary  of  the  domain  of  integration. 
Efficiency  was  achieved  by  the  judicious  application  of  the  alternating  sign 
property.  In  Table  VII,  53  new  multidimensional  quadrature  formulas  with  boundary 
partial  derivative  correction  terms  were  given.  Numerical  results  for  double  and 
triple  integrals  indicate  the  new  formulas  are  both  accurate  and  efficient. 

In  applications  where  first-  and  perhaps  mixed  second-order  partial  derivatives 
of  the  integrand  are  readily  computed,  the  new  formulas  should  be  preferred  to 
traditional  rules.  Specifically,  some  of  the  new  formulas  may  be  useful  in  finite 
element  applications. 

Future  studies  should  consider  the  addition  of  weight  functions,  more  general 
domains  of  integration,  higher-order  partial  derivative  correction  terms,  additional 
and  arbitrary  node  locations,  and  the  use  of  Gregory-type  difference  correction 
terms  in  place  of  the  first-order  partial  derivative  correction  terms. 

Finally,  we  conjecture  that  formula  EM  143  (Table  VTI)  may  be  the  first  of  a 
class  of  Gaussian-type  partial  derivative  corrected  multidimensional  quadrature 
formulas. 
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